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Abstract 

We propose a discontinuous finite element method for small strain elasticity allowing 
for cohesive zone modeling. The method yields a seamless transition between the 
discontinuous Galerkin method and classical cohesive zone modeling. Some relevant 
numerical examples are presented. 


1 Introduction 


In this paper we develop a discontinuous finite element method for cohesive 
zone modeling using the approach hrst suggested by Hansbo and Hansbo [3] . 
Unlike in the standard pre-failure treatment of cohesive zones, which consists 
of tying the meshes together using a penalty approach, we use a combination 
of Nitsche’s method and the cohesive law governing the interelement stiffness, 
thus allwoing the same discretization method in both pre-failure and post¬ 
failure regimes. This means that the method is consistent with the original 
differential equation and no large penalty parameters are required for accurate 
solutions even in the pre-failure regime. The approach was implemented for 
cohesive cracks by Heintz and Hansbo [1] in an XFEM setting, but here we 
consider a discontinuous Galerkin method allowing for discontinuities appear¬ 
ing only between elements. 

An approach similar to ours has been suggested by Mergheim, Kuhl, and Stein- 
mann [7] , an later used by Pretchel et ah [S] and Wu et ah m- The method of 
[7] however uses a different blending of Nitsche’s method and cohesive zones. 
There a discontinuous Galerkin method is used only in the pre-failure regime 
and a switch to a standard cohesive zone approximation is performed at a 
given traction threshold. To ensure a continuous transition between the dis¬ 
cretization methods, a matching of discrete tractions between the two cases is 
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performed. This matching is cumbersome in a more general situation of non¬ 
matching meshes across the cohesive zone. In this paper we avoid this switch 
and a more generally applicable method results. 

An outline of the remainder of the paper is as follows. In Section 2 we dehne our 
discrete method in a linear setting; in Section 3 we discuss and motivate the 
cohesive law that we favour and the resulting secant compliance we use in our 
numerical solution process; and in Section 4 we give some numerical examples 
of our approach. Finally, in Section 5, we give some concluding remarks. 


2 The model problem and discretization method 


2.1 Linear elasticity with a single cohesive interface 


W consider first an elasticity problem in Usd = 2 or 3 dimensions with a smooth 
boundary F dividing hi into two parts hli and 122- The displacements = [«*]£! 
has restrictions to the different domains Sj = u\q^, and we denote by |s] = 
Si|r — S 2 |r and let n denote the outward normal vector to dQ and to hli on 
F. Then, a linear elasticity problem with cohesive layer F can be written: Find 
u and and the symmetric stress tensor cr = such that 

cr = A V ■ ul + 2^e{u) in f2i U 122, 

— V ■ cr = f in hli U 122, 
u = q on (9r2D, 

( 1 ) 

(T ■ n = h on diljq 
|cr ■ n] = 0 on F 
|s] = —Kcr ■ n on F 


Here A and p are positive constants called the Lame constants, satisfying 
0 < pi < /i < p 2 and 0 < A < cx), and e (s) = is the strain tensor 

with components 


1 / dti,- (9 m 1 




u) = - 
^ ’ 2 


+ 


dxi dx. 


Furthermore, V cr = 
Sij = 0ifiy^j,f and 


Ej=i daij/dxj I = with 5ij = liii = j and 

h are given loads, g is a given boundary displacement, 
and n is the outward unit normal to dQ. Finally, K is a symmetric positive 
semi-dehnite flexibility matrix (constitutive law on F). For example, with 
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isotropic elasticity on F we have that 


K = a I + {/3 — a)n ® n, or Kij = a6ij + {(3 — a)ninj, 

where ® denotes outer product, with a > 0 and (3 > 0 denoting the compli- 
cancy in the direction tangential and normal to F, respectively, cf. [3j. (In this 
paper, a more general compliance, with cross coupling between normal and 
tangential directions will be considered.) 


2.2 A discontinuous Galerkin method for linear eohesive zones 


Consider a subdivision of hi into a geometrically conforming hnite element 
partitioning = {T} of kl. Let 

P^{T) = {?;: each component of is a polynomial of degree < fc on T}, 


= {ve ; v\t e [p^{T)f^^ VT e T'^}. 

We also introduce the set of element faces in the mesh, (5^ = {F}, and we split 
5 into three disjoint subsets 


5^ = 5/ U U Jat, 

where is the set of faces in the interior of hi and and are the sets 
of faces on the Dirichlet and Neumann part of the boundary, respectively. 
Further, with each face we associate a hxed unit normal n such that for faces 
on the boundary n is the exterior unit normal. We denote the jump of a 
function v G at an internal face F G '^j hy |t;] = v~^ — v~, and |t;] = v~^ 
for F G ^D, and the average (v) = (v+ + v~)/2 for F G and {v) = for 
F ^ do, where = lim^on(a; Ten) with x E F. 

For the modelling of cohesive interfaces, we here assume that the solution may 
be discontinuous across each element face F, and thus the role of F in ([^ is 
now taken by all element faces. 

The DG method can then be formulated as follows: Seek G such that 
Ohiu^,v) = Lhiv) for all r; G (2) 
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The bilinear form is defined by 


ah{u^,v)= ^ £{v)dx 

Ter'* ^ 

- Y. [ {(^(u^)-n) ■ {lvj+K {(t{v)- n)) ds 

- Y f {(t{v) ■ n) ■ + K ■ n)) ds ( 3 ) 

F&SiUSd ^ 

+ Y [ ■ '^) ■ (-^ (^('“^) ■ '^)) 

Fe^iUdo 

FeSiUdo 

and the linear functional by 


Lh{v)= f f-vdx+ [ h-vds- 

FeSjv F&So 

+ H / {ShO) ■ {v + K(t{v) ■ n) ds. 

f&^d 


■ n ■ g ds 


(4) 


Here Shis a matrix which depends on the interface conditions of the problem, 
the local meshsize, and a penalty parameter 7 := {2g + 3 A) 7 o, where 70 is 
a dimensionless number which has to be large enough for the method to be 
stable. The stability of the method increases with increasing flexibility, so the 
choice of 70 needed for stability in the case of zero flexibility can be used in all 
other cases (numerical values for 70 can be found, e.g., in 0). More precisely, 
on a face F with diameter hp, 

Sh\F=(^I + K) (5) 


On each face F, the mesh parameter is dehned by 

^ f (^meas(T+) + meas(T“)^ /2meas(F) for F C 9T+ fl dT~, 
[meas(T)/meas(F) for F C 9T fl 


We note that as the flexibility goes to zero, we approach a standard discon¬ 
tinuous Galerkin method for elasticity. Looking instead at the limit case of 
h —0 (assuming K is invertible) we retrieve a standard formulation for co¬ 
hesive laws where the only term contributing to the stiffness matrix from the 
interfaces is the interface stiffness term 


E 

FeS/USi) 



(k-‘[u'‘]).H ds. 


4 


The proposed method thus seamlessly blends discontinuous Galerkin with 
standard FEM for cohesive interfaces. 

By use of Green’s formula, we readily establish that the method (|^ is consis¬ 
tent in the sense that 

ah{u-u^,v) = 0 (7) 

for all V G and for sufficiently regular, which is key to retrieving 
optimal accuracy of the method. Stability follows from the analysis in [2ll8ll5] . 
We also mention the work of Juntunen and Stenberg [6], where an analysis of 
this approach for handling general boundary conditions for Poisson’s equation 
is given. 


3 Cohesive law implementation 


We are now interested in the case when the interface compliance depends 
on the jump of the solution, K = FCdu]). In the numerical solution of the 
cohesive zone FE model, we replace this compliance by a corresponding secant 
compliance as follows. 

The compliance between the elements can be dehned a priori by cohesive zone 
models. Frequently, cohesive zone models that are easy to implement are cho¬ 
sen to model the initiation of cracks. These models are often un-coupled, mean¬ 
ing that there is no relationship between the normal and tangential stresses 
other than the fracture criterion. However, it is reasonable to imagine that the 
tangetial stiffness is effected by the reduction of the normal stiffness. Thus, 
in order to couple the cohesive behavior of the interfaces in mixed mode, we 
choose to derive the traction-separation laws from a weighted energy release 
rate surface, cf. m- 

We denote the energy release rates in pure normal and pure tangetial direc¬ 
tions Fi {un) and Tj: (ui), respectively, where, for convenience, the normal and 
tangetial jumps in displacement are denoted by Un = |w] ■ n and Ut = Iw] ■ t, 
where t is the tangent vector to the given face, such that n and t constitute 
a right-handed ON system. The energy release rates are obtained from the 
interface traction on according to 


‘^n 

Fi (un) = Fi {un, 0) = J (T («„, 0) ■ ti dw^, (8) 

0 

Ut 

Fii {ut) = Fi (0, Ut) = j cr (0, Ut) ■ t dut- (9) 

0 
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By use of a polar coordinate system, a dimensionless effective separation A 
can be defined together with an angle ip that determines the mode mix. The 
mode mix and the effective separation are dehned as 


ip = arctan 


/ UncUn 
\ U^pUt 


( 10 ) 



Here, u^c and Utc are the critical normal and tangential separations in pure 
modes. The normal and tangential separations, Wn/wnc and Ut/utc , are dehned 
as the projections of the effective separation on each respective pure mode axis, 
cf. Fig. 



Fig. 1. Illustration of effective separation and mode-mixity 
It then follows that Un and Ut are given by 


Un = XUnc COS (ip) , 


( 12 ) 


Mt = Autc sin ((p). (13) 

In order to obtain the complete contributions of the energy release rates in 
each pure mode and not the projections, some additional dehnitions. Tin and 
Tit are introduced. For example, we may choose to dehne Tin ■= Aunc and Tit '■= 
Xutc- As a hrst step in the development of the cohesive law, two independent 
functions are htted to experimentally measured energy release rate curves, see 
Fj and Fu in Figj^ 
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Fig. 2. Schematic illustration of the fitted energy release rate curves in normal and 
tangential directions. 

The shapes of the traction-separation curves in each pure mode, respectively, 
are obtained by differentiating the energy release rates in each pure mode with 
respect to each pure mode relative separation, Un and Ut. From these curves, 
laws are chosen that captures the most essential features of the curves. Figure 
[ 3 ] shows two idealized schematic curves. 

In order to capture the behavior of the cohesive law in mixed mode, the two 
energy release rate curves in Fig. [^are combined to yield a surface where the 
axes are total energy release rate, F, relative normal and relative tangential 
separations, and Ut respectively, see Fig. 



Fig. 3. Schematic illustration of normal (solid) and tangential (dashed) traction-sep¬ 
aration curves. 

The surface representing the weighted energy release rate, F (A, (p) is generated 
by a weighted sum of the experimentally determined energy release rates in 
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pure normal, il, and pure tangential, /ji, directions according to 


r, (A, ^) = f {,p) A (A) ul + (l-f (f)) A, (A) ul (14) 


where / (cp) is the weight function. 

The stresses for any given mode mix are given by partial differentiation of F 
with respect to each relative separation, Un and Ut, respectively. 


dr _ dr dX dr dip 
dUn dX dUn dip dUn ’ 


(15) 


dr _ dr dX dr dip 
dut dX dut dip dut 


(16) 



V (mm) (mm) 


Fig. 4. Weighted potential surface. 

The secant compliance is computed as follows. We hrst establish the secant 
stiffness matrix St as 

(T ■ n cr ■ t (j ■ t cr ■ n 

St = - n <8) n H- 1 ® t H- n ® t H- 1 ® n, 

Un Uf Un Uf 

followed by computing the secant compliance as Kt = S^^. The interface 
stiffness. Shir, is then given by 


Ship = 


hi 


I + Kj 


-1 


7 


(17) 




























4 Numerical example 


A specimen with two inclusions and an initial crack, see Fig. is used as 
a simple example to show the applicability of the modeling technique. The 
dimensions of the specimen are given by; W = H = 1.00 mm, D = 0.20 
mm, a = 0.20 mm. The lower right inclusion is located at center coordinate 
(0.75,1.00) mm and the top left inclusion is located at (0.45,1.10) mm. The 
crack is located at center coordinates (0.40,0.90) mm and it is inclined at an 
angle of 33° to the horizontal axis. The boundary conditions for the specimen 
are set to be clamped on the bottom edge, i.e. Ux{x,Q) = My(a:, 0) = 0. The 
top boundary is constrained horizontally Ux{x, 2H) = 0 and the displacement 
is controlled vertically Uy{x, 2H) = A, see Fig. Two different set-ups are 
modeled for comparison. The hrst is a specimen where the inclusions have 
the same material properties as the rest of the specimen with elastic material 
properties; B = 10 MPa and z/ = 0.45. The second is a specimen where the 
Young’s modulus of the inclusions is 100 times greater than in the rest of 
the specimen. The maximum cohesive strengths are set to 1 MPa and the 
maximum critical separations are set to 0.02 mm in the cohesive sawtooth 
model giving a fracture energy of 0.1 J/mm^. Note that these properties are 
the same for both set-ups. 

One of the major issues with this type of modeling is mesh dependency. How¬ 
ever, if a large number of elements is used the mesh dependency is obviously 
reduced. Furthermore, the compliance between all continuum elements in¬ 
troduce numerical issues which can be reduced by an increase of the elastic 
stiffness of the cohesive zone model sufficiently to minimize the compliance. 



Fig. 5. Dimensions of the Single Edge Notched specimen. 

In the present model, the compliance is allowed to be initially zero and then 
gradually increase as the load is increased. Damage initiation, and essentially 
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crack propagation, is enabled by a decrease of the stiffness according to (17) 
where the interfaces, as stated in the dehnition of the method, are given as 
the boundaries between all the continuum elements (this is of course not a re¬ 
quirement, as a mix of continuous and discontinuous methods is also possible). 
Thus, cracks are free to form, nucleate and propagate along the continuum 
element boundaries by Nitsche’s method instead of the standard approach of 
using cohesive elements, as in, e.g., m. 


For the hrst set-up, see Fig. the crack initiates as expected and then it 
propagates without considering the inclusions. In the second set-up, however, 
the crack is arrested by the stiffer inclusion boundary and deflects downwards 
around the lower right inclusion to hnally to continue to propagate to the free 
edge of the specimen. It can be seen for both set-ups that there is virtually no 
compliance issues prior to any cracks forming. Both simulations, thus shows 
the applicability of the modeling technique. 



Fig. 6. Increasing deformation from left to right for the first set-up. 





] 




0. 

1 



i 



Fig. 7. Increasing deformation from left to right for the second set-up. 
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5 Concluding remarks 


In this paper, we have suggested an FE method which seamlessly blends the 
discontinuous Galerkin method with classical cohesive zone models. There is 
no need for interface elements as the interelement stiffness is represented by a 
modihcation of the weak form. There is no need to identify threshold values 
for transitions between discretization approaches since the same bilinear form 
is used for all cases of interface stiffness. The method also directly allows 
for modeling cohesive zones between non-matching meshes, unlike the similar 
approach suggested previously in [7], which does not immediately generalize 
to this case. 
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